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Abstract 

The new concept of numerical smoothness is applied to the RKDG (Runge-Kutta/Discontinuous 
Galerkin) methods for scalar nonlinear conservations laws. The main result is an a posteriori error 
estimate for the RKDG methods of arbitrary order in space and time, with optimal convergence rate. 
In this paper, the case of smooth solutions is the focus point. However, the error analysis framework is 
prepared to deal with discontinuous solutions in the future. 

1 Introduction 

The Runge-Kutta/Discontinuous Galerkin (RKDG) methods are among the most popular modern numer- 
ical methods for nonlinear conservation laws. Due to the complexity of the schemes and the nonlinearity 
of the problems, error analysis theory for the RKDG methods is not yet satisfactorily completed. A brief 
summary of the currently available error analysis results can be found in the recent papers by Q. Zhang 
and C.-W. Shu, [9] and [lOj. Here in this paper, a new error analysis is given based on the innovative 
concept of numerical smoothness. The main result of this paper is a practical a posteriori error estimate 
of optimal order , which depends on a set of computed smoothness indicators. 

In most a priori error analysis of time-dependent problems, local error is referred to original PDE 
solutions to take advantage of their smoothness, consequently error propagation has to be referred to 
numerical schemes (e.g. [TO]). In a posteriori error analysis, local error can only be referred to numerical 
solutions, consequently error propagation is referred to PDEs. Since numerical solutions of PDEs are 
typically not smooth functions (discrete point values in finite difference schemes, piecewise polynomials in 
finite elements, etc.), when local error is referred to numerical solutions, it is often given as residuals, such 
as in the well-known duality method [T] [2]. In our a posteriori error analysis of the scalar conservation 
laws and RKDG, we use the Li-contraction between the PDEs' entropy solutions for error propagation 
analysis, and rely on numerical smoothness instead of residuals to estimate local error. 

The idea of using numerical smoothness in the error analysis of nonlinear hyperbolic conservation 
laws is a migration of the idea of using numerical smoothing in the error analysis of nonlinear parabolic 
equations solved with complex schemes [?][§]■ For nonlinear equations solved with complex schemes, 
we base the concept of numerical smoothness on a set of efficiently computable smoothness indicators. 
When the indicators remain bounded during actual computation, we consider the numerical solution as 
being numerically smooth. Due to the equations' nonlinearity and the schemes' complexity, we usually 
cannot give an a priori proof the boundedness of our smoothness indicators. However, we can always 
compute the indicators along a numerical solution. We define the indicators at t n in such a way that 
we can prove a local error estimate of optimal rate for the time step [t n ,t n +i], where the indicators play 
the role of high order derivatives as in most a priori error estimates. It will be shown that numerical 
smoothness indicators deliver much more abundant information then residuals. Consequently, we can 
get better error estimates and more useful information toward adaptive algorithms. 

Numerical solutions are not always numerically smooth. The usual measures taken for the purpose of 
achieving numerical stability are actually also what is needed to achieve numerical smoothness, because 
one of the main ingredients of these measures is numerical diffusion. Smoothness indicators serve two 
purposes: (1) to watch the smoothing and/or smoothness maintenance performance of the scheme; and 
(2) to provide smoothness information for local error estimates. For the RKDG schemes, we use the 
Godunov upwind flux, the TVD-RK schemes and a strengthened CFL condition. After taking all these 
measures, it is extremely hard to prove the boundedness of our smoothness indicators. However, the 
importance of the whole idea resides on the fact that we can use the computed smoothness indicators 



to circumvent the difficult proof of numerical smoothness, and move on to prove sharp error estimates. 
For complex nonlinear problems, this kind of circumvention is probable the only way to achieve practical 
error estimates. 

While doing the proofs, we realized another advantage of the numerical smoothness approach. Since 
we are working on DG finite element solutions and the entropy solutions just evolving away from those 
finite element solutions, we have easy access to the Loo-norm estimates, which in turn gives us L\ and 
Z/2 estimates. For error propagation, L\ contraction is the best tool. For the finite element formulations, 
Z/2-norm is natural. For the nonlinearity, Loo-estimates are crucial. Having access to all three, we are 
able to do the error analysis of the finite element methods of the nonlinear problems, where the global 
error estimates do not have an exponentially growing factor. 

In the solutions of nonlinear conservation laws, there are shocks and contact discontinuities. In the 
discontinuous Galerkin finite elements, there are the technical discontinuities of the piecewise polynomials. 
The analysis of this paper is limited to the case of smooth PDE solutions; henceforth, only the technical 
discontinuities are treated. While we do not assume anything directly on the smoothness of the PDE's 
solution, we only consider the case that all the components of our smoothness indicators are well-bounded. 
In fact, the boundedness of the smoothness indicators indicates that a smooth PDE solution is being 
approximated. Our smoothness indicators are capable of detecting shocks and contact discontinuities 
(including high order discontinuities) [3]. Our Li-contraction error propagation analysis remains valid 
in dealing with shocks and other discontinuous solutions. However, we restrict ourselves to the case of 
smooth solutions in this paper. Many discontinuous solutions of conservation laws are piecewise smooth. 
The work of this paper can also be considered as analyzing error in a smooth piece of a discontinuous 
solution. Clearly, we need to understand how to obtain optimal error estimates on smooth pieces of 
solutions, before we focus on the error analysis at shocks and contact discontinuities. In this sense, this 
paper is the first step of the project of analyzing the error of RKDG methods with numerical smoothness. 

In this paper, our goal is to show the new error analysis ideas and the nature of the results. In order 
to focus on the framework, we do not trace all the constants involved in the error estimates. Instead, we 
show how they should be computed with enough details to reveal their dependency, computability and 
boundedness. A separate technical report will be prepared to show the fine details. Since some of the 
constants depend on the flux function /, certain details are better shown with numerical examples. No 
generic constant will appear in this paper. 

The nonlinear conservation law problems and the RKDG schemes are well-known. For a survey 
article, see the lecture notes [5] by C.-W. Shu. Consider the one-dimensional nonlinear conservation law 

«* + /(«)« = (1) 

in a bounded interval f2 = [a, b]. In order to focus on the new ideas and the new tools of the proof, we 
stay with the simple case of west wind ( f'(u) > ) . Let the initial condition be 

u(0,x)=ui(x) (2) 

and the upwind boundary condition be 

u(t, a) = u L (t). (3) 

Assume that the flux function f(u) is sufficiently smooth and the initial and boundary conditions are 
smooth and consistent to guarantee that the entropy solution u(t, x) is smooth near x — a for all t > 0. 

Partition Q with a — X-1/2 < X\n < • • • < x m -l/2 = b. Let h — Xj+1/2 — Xj_ 1 / 2 be the same for all 
cells flj — [xj-1/2, Xj + i/ 2 ]- To solve the problem with the discontinuous Galerkin method, we take the 
standard discontinuous piecewise polynomials space Vh- When the degree of a local polynomial is up to 
p, Vh ~ {v € L 2 {Vi) : v\n - G n p }, where n p is the set of all the polynomials of degree less than or equal 
to p. In each cell flj, a semi-discrete solution uh satisfies 

{u h , t ,v) nj = (f(u h ),v x ) nj + f(Uh(xJ_ 1/2 ))v(xf_ 1/2 ) - f(uh(xJ +1/2 ))v{xJ +1/2 ). (4) 

Here the Godunov flux is employed under the west wind assumption (for simplicity). At the upwind 
boundary X-1/2 = a, we set 

u h (t,X-i/ 2 ) = u(t,a) = Uh(t). (5) 
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At the initial time t = 0, Uh(0,x) is taken to be the Z/2-projection of ui(x). 

For temporal discretization, we take a standard TVD-RK scheme of order k [5]. For example, when 
k = 3, in the time step [£ n ,£n+i]> we compute the fully-discrete solution u n+1 £ Vh from u n by the 
following scheme. With the notation 

Kj{u,v) = (f{u),v x ) nj + f(u(xJ_ 1/a ))v(xf_ 1/2 ) - f(u(x~ +1/2 ))v(x- +1/2 ), 

the scheme is that, for r n = t n +i — t n and any v G Vh, 

(u„ 1 ,v)n :l = (u n ,v)n ] + T n Hj{u n ,v), (6) 

{u c f,v)n j = t(i4>«)oj + tCmS.' 1 ,^)^ + '^-'Hj(u^ 1 ,v), (7) 

« +1 ,ii) f2j . = ^(u n ,v)n 3 + ^{u^ 2 ,v)n 3 + 2 ^ / H j {u c r ; 2 ,v). (8) 

At t = 0, we make u§ = u^(0, s). The upwind boundary values of it^' 1 , u^' 2 , and w n+ i are taken according 
to u L {t). 



2 Error propagation and numerical smoothness 

The PDE's entropy solution satisfying the original initial condition is denoted by u(t,x). Throughout 
this paper, we only consider one numerical solution, namely the computed numerical solution, which is 
denoted by u^j as above. In order to present the local error analysis in a time step [t n ,t n +l], we use 
the semi-discrete solution that passes (t n ,M n ). For the briefness of notations, we use u h (t,x) for this 
temporally piecewise semi-discrete solution, which has a new initial value in each time step. u h takes the 
upwind boundary value given in |5|. Since we do not simultaneously work on the local error analysis of 
two different time steps, the notation u h (t, x) should not cause ambiguity. In each time step, we also need 
the PDE's entropy solution which passes (t n ,u^). We denote this entropy solution by u(t,x). u(t,x) 
satisfies the upwind boundary condition |3|. Of course, u is also defined piecewise in time. The following 
error splitting diagram may help the reader in remembering the notations for these solutions. 







u(t„+i) 






u(t n+ i) 






u h {t n+1 ) 




u(t„) 















In the diagram and also in the rest of the paper, sometimes we hide one of the two independent 
variables in the notation of a solution to make the expressions shorter. 

The error analysis of this paper is based on the error splitting in the diagram. In order to estimate 
the global error u(t n+i) u n _^_i at time we split it into three parts as shown in the diagram. 

\\u(tn+l) — U n +l\\ < ||w(tn+l) — "(tn+OH 
+ ||«(t»+l) — u h (t„ + i)\\ 

+ \\u h (t n+ i) - (9) 
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The first part u(t„+i) — u(t n +i) is the propagation of the global error u(t„) — u n by the PDE. Due 
to the Li-contraction property of the scalar conservation laws, since u and u satisfy the same upwind 
boundary condition (J3j) , we have 

||u(t n +i) - u(t n+1 )\\ Ll{n) < \\u(t n ) - u n \\ Ll{sl) . (10) 

The second part of the split error is the local spatial discretization error u(t n +i) — u h (t„+i). Since u n 
lives in a discontinuous finite element space, u is certainly not smooth in the classical sense. In fact, the 
discontinuity of u(t n ) = u„ at Xj_i/ 2 will travel into the cell fi,-. Thus, at any time t £ (t n , t n +i], there is 
either a shock, a contact discontinuity, or a rarefaction wave of u in Qj. However, if the solution u(t, x) 
is smooth around Qj, we intuitively know that the discontinuity of u is only technical and it must be 
very tiny, and u must be smooth away from its discontinuities. We will quantitatively substantiate the 
intuition in the definition of the spatial smoothness indicator. Then we will use the indicator to estimate 
the spatial local error. 

The third part of the split error is the local temporal discretization error u 1 (i n+i) ^n+i* To estimate 
this part of the error, we will need the temporal smoothness of u h for t £ \t„, t n +i\. Since u h is an ODE 
solution with initial value u„, we can also establish the needed smoothness. 

Here it is to be noticed that the analysis relies on the smoothness properties of u and u . Since 
both of them have u n as their initial value, the smoothness level of them depends on how u n has been 
computed. In other words, some kind of numerical smoothing or smoothness-maintenance should have 
been built in the scheme. In this paper, we do not intend to prove such smoothing or smoothness- 
maintenance ability for the RKDG methods. Fortunately, the computed smoothness indicators in our 
numerical experiments show that the RKDG methods do have the desired ability to keep a numerical 
solution "smooth" (when/where the solution should be smooth). We only prove error estimates by using 
the smoothness indicators. 

3 The smoothness indicators 

In order to rigorously and quantitatively define the concept of numerical smoothness for the RKDG 
method, we define the following spatial and temporal smoothness indicators for each time step [t n ,t n +i]. 

• Spatial smoothness indicator: S% = Sp(u^), 

• Temporal smoothness indicator: T„ — Tk(u n ). 

Here S stands for space, T stands for time, p is the degree of the polynomials in each cell, and k is the 
order of the Runge-Kutta scheme. 

3.1 Definition of T% 

The temporal smoothness indicator consists of the temporal derivatives of u h at t = t n . Namely, 

k c h h d k+1 u h 

The first derivative ui(t„) is computed as in the implementation of the forward Euler scheme 

(u t (t n ),v)n j = Uj(u c n ,v). 

The formula for computing the second derivative can be obtained by taking derivatives with respect to 
time on both sides of the semi-discrete scheme Q: 

= (f'(u h )ut,v x )n j 

+ f'(u h (xJ_ 1/2 ))ut(x~_ 1/2 )v(x+_ 1/2 ) 

- f'( uh ( x J+l/2)) U t( x J+l/ 2 ) V ( x J+l/2)- 

To compute Utt(t n ) with this formula, on the right hand side, we replace u h by u n and Ut by the computed 
first derivative Ut(t n ). The high order derivatives can be computed similarly. 
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The ability of the indicator to reveal numerical solutions' smoothness, discontinuities, and possible 
numerical "instability" phenomena has been reported in 4.. Since contains the initial temporal 
derivatives of u , it can be used for the temporal local error estimation without any transformation. 



3.2 Definition of S* 

The spatial smoothness indicator _JJ contains not only the spatial derivatives of M^j within each cell, but 
also the jumps of the derivatives across the cell boundaries. Namely, for the cell flj = [xj-i/2, 2^+1/2], 

Sl tj = {M° ntj ,Ml^--- ,M%j,D°j,Dn t j,--- 

where 



and 



1 l c 4 1 d l 



J l — l\,-f l — T l — D ! JjP+i+M-ifl+a) 



for some properly determined constants a £ [0, 1) and /i £ [0, 1]. 

When j — 0, L l n needs to be defined separately. To this end, we use the upwind boundary function 
u(t,x^ 1 / 2 ) = UL(t). By calculation from the conservation law |TJ), we must have 

L° nfi = u(t n , x_ 1/2 ) = u L (t n ), 

Ln -°- Ux{tn ' x - 1/2) -~f(u L (t n )y 

T 2 _ ,. s_ 2f"(u L (t n ))[f t u L (t n )] 2 - f'( UL (t n ))^U L (t n ) 

L n>0 - u xx (t n ,x_ l/2 ) - _____ , 

and so on. For later use, we extend ujj to f2_i = [__i/2 — h,x_ 1 / 2 ] by u£ — u(t n ,x), where u(t n ,x) is 
obtained by a short time tracing back from u(t,a) = Under a proper smoothness assumption on 

UL(t), such tracing back is well defined. By Taylor expansion, 

Un[x) = u(t n , x) = L% + Ln >0 (x - k_i/ 2 ) H h L^ (x - x_ 1/2 ) p /pl + R n ,o(x), x e Cl-l. 

The residual R n ,o(x) = 0((x — x_ 1 / 2 ) p+1 ) is of higher order. Given a smooth boundary function ul, one 
can determine a constant D, such that 

\R n , (x)\<D\x-X- 1/2 \ p+1 /(p+l)\ 

In other cells (j > 0), let R n ,j(x) = 0. The expansion part of is not computable, it only lives in the 
proof. The expansion is defined in this way to be consistent with the boundary condition satisfied by u. 

It is obvious that the values of M l n j and L n j should be of 0(1), unless there is a shock or contact 
discontinuity somewhere around Q.j. It is also easy to guess that the jumps J n j should be small, 
otherwise the numerical solution may have lost too much smoothness around the cell boundary. How 
small should the jumps J l n j be? Both our error analysis and numerical experiments suggest that D l n j = 
Jnj //i p+1+M_ '' 1+a! ' should be at most of 0(1), unless there is a shock or high order discontinuity within 
or near the cell. This is the reason for having D l „ j instead of J l n j serving as a part of the smoothness 
indicator. 

How is a determined? It is well known that, with high degree DG elements, the time step size r n 
should satisfy a strengthened CFL condition of the form 

.l+a 
T n < jh 

In [9], for example, a = 1/3. For the definition of D n j, I = 0, ■• ■ ,p, we need a = /i/p. In fact, since 
is the jump of the piecewise constant function to match the total variation of -§*^u, the average 

of the jumps J p , must be of 0(h). That is, we need p + 1 + /i — p(l + a) = 1 in the definition of D v n ■, 
or equivalently fi — pa, to have = 0(1). 
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The amount of work for computing SJJ is proper. In fact, S£ contains the minimal amount of smooth- 
ness information for us to estimate the optimal approximation error of the piecewise polynomials of degree 
p to the PDE solution. The amount of work for computing T„ is also proper, for the same reason. It 
might be possible to estimate from S£ (if k and p are related in certain way), but a directly computed 
T„ should sharpen the temporal local error estimate. 

We consider the numerical solution as a good approximation of a smooth PDE solution if and only 
if S„ is reasonably bounded. It will be a future issue to study how to classify the numerical solution in 
case S„ is not reasonably bounded. We will have to distinguish different patterns of SJJ. What indicates 
a well-caught shock, or well-approximated transition to a shock? What indicates a well-approximated 
high order contact discontinuity? What indicates numerical "instability"? How to adaptively deal with 
each of these cases? The temporal smoothness indicator T% should also be studied for the same issues. 



4 The main error estimates 

Theorem 4.1 Let u(t,x) be the entropy solution of the nonlinear conservation law satisfying the 
initial condition and upwind boundary condition J5|). Let be the numerical solution computed by a 
TVD-RK-DG scheme with piecewise polynomials of degree p and the TVD-RK scheme of order k, on the 
partition of Q, described in Section 1. Assume that u and are bounded by a constant U in [0,T] x f2. 
Let P = maxi w i<(7 f'(w). Assume that the time step size r n for each step satisfies the standard CFL 
condition Pr n < h and the strengthened CFL condition r n < r yh 1+a , for a constant fi £ [0,1], a positive 
constant 7, and a = p/p. 

Lf there is a positive real number M, such that, for all t n < T, all the components of S% and are 
bounded by M, then the spatial and temporal local error in [t„,t„+i] satisfy 

\\u(t n+ i) - u h (t n+1 )\\ Ll(n) < r n h p+ ^(S V n), (11) 

||« (tn+l) - «n+l|Ul(n) < Tn + 1 Q{Tn, S%) , (12) 

where J-(S%) and Q(T„, S„) are computable functions of the indicators. As a consequence of the error 
splitting |5|), the L\- contraction property \lff\j , and the local error estimates flZ]) and [!§), 

\\u(t n +i) - u c n+1 \\ Ll(n) < \\u(t n ) ~ u c n \\ Ll(ft) + r n [/i p+M J r (5'^) + T n G(T%, £%,)]. 
Finally, at the end of the computation (tn = T), 

\\u{T)-u c N \\ Ll(n) < HO) -ug|U l(n) + ^ 1 r n [h p+ ^(Sl) +r*g{TtSl)]. (13) 
Proof. It suffices to prove and \12\ . The next two subsections will carry the proofs of (jTTj) and 



( 12 1 respectively. # 



Remark. In the literature, /i — 1 is considered to be the optimal convergence rate. We keep p, as a 
parameter to cover those possible non-optimal cases. However, when the initial solution is smooth, we 
always have y, = 1. For p > 3, a = ji/p is not too restrictive. There is no restriction on 7 in the proof, 
although 7 will appear in the function J-^S^). The actual restriction on t„ is in real computation. If r n is 
too large, the RKDG scheme fails on numerical smoothness maintenance. See the numerical experiments 
in Section 5. 



4.1 Estimating u(t n+ i) —u h (t n+ i), proof of (11) 



We begin with introducing an auxiliary piecewise PDE solution u e . First define a local strong solution 
Uj of the conservation law. The initial values of Uj are given on the line segment {£„} x (flj-i U fij) by 

Uj(t„, x) = M°j + Mnj[x - ij-i/a) H h M% tj (x - x 3 _ l/2 Y /p\ j = 0, 1,- •■ ,ra - 1. 

It is easy to see that, in Qj, Uj(t n ) = u c n \ in Qj-i, 

Uj(t n ) = Un + J°,j + Jn,j(x ~ Xj-l/2) H + Jn,j(. X - x J-l/2) f '/p! ~ R»,j(x). 
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As a strong solution of the Cauche problem of the original conservation law (Tl|), Uj certainly exists in the 
region TZ„j = {(t, x)\t g [t n , t n +i],x G Qj-i U Qj,x < x j+1 / 2 ,x — x + f'(Uj{t n , x))(t — t n )}. This is the 
trapozoidal region covered by the characteristic lines originating from Clj-i U Qj. When t„ satisfies the 
standard CFL condition f3r n < h, it is easy to see that [i n ,in+i] x ®,j C TZn,j C [tnjtn+i] x (^j-i Uf! 3 ). 

At the upwind boundary, let U-i = u for x G = [x-1/2 — h,X-i/2\- Due to the smoothness of 
UL(t), one can determine the value of u in f2_i by tracing back (but not computable). Now, we are ready 
to define the local piecewise PDE solution by 

u e (t,x) — u e j{t, x), (t,x) G [t n ,t n+1 ] x Clj, j = -1, 0, 1, • • • ,m - 1. 

Since Uj(t n ) is a polynomial in fij-i U Qj, u e is smooth in [t n ,t n +i] x Qj for sufficiently small r n . To 
reveal more details on the smoothness of u e , we have the following Lemma. 

Lemma 4.2 There are constants N„ j (I = 0, 1, • •■ ,p + l,j = 0, 1, •• ■ , m— 1,), which depend on the flux 
function f and can be computed from M®j, M*j, • • • , M p j, such that, 

gl 

\\- d — i u e j {t,x)\\ tooi T Ln ^<N nij , 1 = 0, 1,--- ,p, j = 0,1, ••• ,m-l. 

Moreover, 

QP+ 1 

II q xP +i u K tn + r ' x ) Hwn f ) < rAr nj\ i = 0,1,- •■ ,m- 1. 

Proof. For the simplicity of notations, in the proofs of this and the next Lemma, we denote the solution 
Uj by w, and by From 

w t + f(w) x = 0, 

we get 

+ f'(w)w^ + f"(w)(w W f = 0, 
™ t (2) + /Wi 2) + 3/"(w)to (2) ™ (1) + f"'(w)(w (1) ) 3 = 0, 

and so on. The boundedness of w is obvious. 

Along each characteristic line, ur 1 ' will not have a blow-up in a short time r G [0, r n ], where r = t—t n . 
The initial value of i^ 1 ' is £tif (tn) = M 1 J +Ml j (x-x^ 1/2 ) + - ■ ■ + M^ j (x-x j _ 1/2 ) p - 1 /(p - 1)!. /"(w) 
can also be computed from Uj(t„). Therefore, u/ 1 ' can be estimated by the values of M° i3 -, M^ j, ■ ■ ■ , M p . 
in the Ltx> norm. Hence iV* j can be computed. 

As for the higher order derivatives, the ODE for each v>"' along each characteristic line is linear in 
and depends on the lower order derivatives w, i// 1 ' , • • • , w^ -1 ' . The initial value of only depends 
on Uj(t n ), therefore the bound N^j of ur ! ' can also be estimated from M°j, M„j, • • • , M p as a result 
of mathematical induction. 

y,(p+i) jjg^g a S p ec i a i property: its initial value is the function (the p+l-st derivative of Uj(t n )). 
Integrating the ODEs about u/ p+1 ^ along each characteristic line, we realize that w {p+1) is proportional 
to r, while the coefficient N p ^ can be computed from , M 1 ^ , • • • , M p ■ . # 

Remarks. Here we make two remarks on the results of Lemma |4.2| (1) Because each ODE is inte- 
grated along a characteristic line for a short time (not longer than a time step), one can simply use M n a 
as a practical estimate for In other words, Nn,j is actually very close to M^j- However, the most 

useful N^j 1 has to be computed through solving the differential inequalities. (2) The factor r in the 
estimate of is crucial for the error analysis later on. It means that, because m| is evolving out of a 

polynomial of degree p, when Uj is approximated by a polynomial of degree p, the error is proportional to 
r. The idea of picking up this r for the local spatial error (in one way or another) comes from reading [10] , 

As it appears in most error analysis of finite element methods, we also need an Z/2-projection of a 
smooth solution. To this end, we consider the cell by cell Z/2-projection of u e into Vh- Denote this 
projection by u p = u p (t,x) g Vh (p stands for projection here), it is given by 

(u p ,v)n j = (u e ,v)n 3 , Vu G Vh- 
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By the Bramble-Hilbert Lemma, the scaling argument, and Lemma |4.2| the following estimates are 
obvious: 

Lemma 4.3 For sufficiently small r = t — t n , 

\\u - u p |k( %) < C,h p+1 \\w ip+1) \\ Ll{nj) < C^rhNlf, 
\W - ^|U 2(%) < C 2 h p+1 \\w^\\ w < C 2 h p+1 rVhN^\ 



and 



K-u p \\ Loo(nj) < C 3 h p+1 \\w^ +1) \\ Lao{(h) < C 3 hF +1 TNlf. 



Here C\,C2 and C3 are the projection error constants in the reference cell. Consequently, in the whole 
domain Q, let N p+1 = maxj , we have 

\\u e -u p \\ Lim <c 1 \n\h p+1 rN p+1 , 

IK - u p \\ L2(n) < C 2 ^W\h p+1 rN p+1 , 

and, for the cell boundary terms, 

(E7 = - V e (x; +1/2 ) - u p (xJ +1/2 )\ 2 )^ < C 3 y/W\h P+1/2 rN p+1 . 

Next we look into the difference u — u e . In the cell Q,j, at time t n +t, both of these entropy solutions 
u and u e depend on their initial value in Qj _ 1 U Qj . More precisely, since /3 = max /' (u) is the maximum 
of wave speed, both entropy solutions u and u e depend on their initial value in [xj-^/2 —fir, Xj-1/2] Uilj. 
Notice that the initial value of these two solutions are the same in flj and the difference of their initial 
values in [aij_i/2 — fi T > x j-i/2] is 

Uj(t n ) - Un = J°,j + J n ,j( x ~ X j-l/2> H r J„,j(x - X j _ 1/2 ) P /p\ - Rn,j{x). 

If D° j, D\ j, ■ ■ ■ , D p n j are bounded by D, r < , yh 1+a , j > 0, and x £ {xj-i/ 2 — /?r, Xj^i/2], then 
K(*n) - <\ ^ \D° n>J h p+1+ ' 1 + J Di J /i p+1+fl - 1 - a ^ 7 /i 1+a + ■ ■ ■ + D p 3 h p+1+ ^- p{1+a) (Pj) p h {1+a)p /p\ I 

< Dh p+1+ » [1 + 07 + • • • + (07) p /pl] 

< De' il h p+1+ ' 1 . (14) 

For j = 0, there is the the extra residual term Rn,o(x), so we have 

|«o(t»)-t£| < Dh p+1+ »{l + l3 1 + --- + (l3j) p /pl]+D(l3 1 h 1+a ) p+1 /(p+l)\ 

< De ''h p+1+ ", (15) 

if Dh (p+1)a ~" < D (which is easy to satisfy). 

Now, by Theorem 16.1 in the textbook [6] by Joel Smoller, we have 

Lemma 4.4 If j3r < h and r < jh 1+a , then 

\\u(t n + r)-u e (t n +T)\\ Ll((lj) < ||<-uf(t n )|| il[x ._ 1/2 _^._ 1/2] < (Pr)(De^h p+1+ »). 
Consequently, 

\\u(t n + r) - u e (t n + T)\\ Ll(m < rh p+ ^De^\n\. 
Remarks. Again, a few remarks may help. (1) Lemma |4.4| takes the inter-cell technical discontinu- 
ities of the numerical solution into account. Obviously, D is playing the role of a smoothness measure- 
ment. (2) Under the strengthened CFL condition, we are able to allow the value of J l n j to be of order 
^p+i+M-z(i+a>) As it is shown in the numerical examples, when the order of the derivative goes up by 
one, the power (ft ) of the jumps goes down by more than one. The strengthened CFL condition seems 
to help here in the error control, at least in the analysis, even if the jumps of the high order derivatives 
grow quickly. When the smoothness deteriorates near the formation of a shock, the strengthened CFL 
condition may play a role of suppressing Runge phenomena, to some extent. Such Runge phenomena 
and their transport to the downstream should be what causes numerical oscillations. 

Now we are ready to state and prove the last lemma to estimate u p — u h , then we will conclude this 
subsection with the main theorem to estimate u{t n +i) — « (tn+i). 
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Lemma 4.5 There is a computable constant Q n , depending on S p , such that 

\\u p (t n + r) - u h (t n + r)\\ Ll(n) < rh p+p Q n . (16) 

Proof. In the cell Qj, u h satisfies the semi-discrete DG scheme that is, 

{ut ,v) nj = (f(u h ),v x )n J + f(u h (x~_ 1/2 ))v(x+_ 1/2 ) - f(u h (x~ +1/2 ))v(x~ +1/2 ). (17) 

Consider the piecewise strong solution tt e , which is the restriction of Uj in Qj. Multiplying u\ + f(u e ) x = 
by a test function v, integrating in Qj, after using integration by parts, we get 

(ut,u)nj = {f{u),v x )n 3 + f(u e (x+_ 1/2 ))v(x+_ 1/2 ) - f(u e {xJ +1/2 ))v(x- +1/2 ). (18) 

Since (u p ,v)n j = (u e ,v)n j for all values of r G [0,r„], (u P ,v)n j = (ut,v)n j . By adding and subtracting 
terms in \1S\ , we get 

(u p t ,v)n 3 = (f(u p ),v x )n 3 +f(u e (x+_ 1/2 ))v(x+_ 1/2 )-f(u p (xJ +1/2 ))v(xJ +1/2 ) 
+ (/(w e ) - f{u p ),v x )a 3 - [f{u{x~ +1/2 )) - ,f(u p (xJ +1/2 ))]v(xJ +1/2 ). 



Now let £ = u p —u h , and let v = f. Subtracting the last equation by (171, we get 

On, = (/(tO-/(Ae„) , 

- [f(u p (xJ +1/2 )) - f(u h (xJ +1/2 ))} £{xJ +1/2 ) 
+ [f{u{x+_ 1/2 )) - f(u h (xJ_ 1/2 ))] e(x+_ 1/2 ) 

+ (/(u e )-/K),f,)o, 

- [/(« <, (x7 +1/a ))-/(u p (x7 +1/a ))]€(a:7 +1/a ). (19) 
First we focus on the third line of the last equation. 



\f(u e (x+_ 1/2 ))-f(u h (xJ_ 1/2 ))\ < 


P\u\x+_ 


1/2) 


hr - 

-U ( Xj _ 


-1/2 


< 


P\u e (x+_ 


1/2) 


e ( — 

- (Sj- 


1/2 


+ 


P \u e (xj_ 


1/2) 


- u p (x~_ 


1/2 


+ 


/3\u p (xJ_ 


1/2) 


hi - 
-U { X] _ 


-1/2 



(20) 

for j > 0. As for j = 0, by verifying u e (a;~ 1 y 2 ) = w (a;^ , 2 ) from the boundary conditions, we have 

\f(u e (xt 1/2 ))- f(u h (xZ 1/2 ))\ < l3\u{xt 1/2 ) -u h (xZ 1/2 )\ 

= p\u*(xt 1/2 )-u*(xZ 1/2 )\. (21) 

Here, we use a brief notation A — u e (x~J_ 1 , 2 ), and B — u e {x~_ 1 ^ 2 ). Recall that u e = vfj in Qj and 
it e = Uj_i in f2j_i. Also recall that Uj is extended to TZ„,j- By the characteristic line theory, we know 
that 

A = Uj(t n + T,Xj-l/ 2 ) = u e j (t n ,Xj_ 1/2 - Tf'{Uj{t n + T,Xj-i/ 3 ))) = u C j(t n ,x J _ 1/2 ~ rf'(A)) 

and 

B = U e j ^ 1 (t n + T,Xj-l/ 2 ) = u e j _ 1 (t n ,Xj_ 1/2 - Tf'(u B j ^ 1 (t n + T,Xj_ 1/2 ))) = u e j _ 1 (t„,Xj^ 1/ 2 - rf'(B)). 
So 

\A-B\ = \uj(t n ,x j - 1 / 2 - rf'(A)) —Uj^tnjXj-i/z —rf'{B))\ 

< \u C j(t n ,x J _ 1/2 -rf'(A)) -Uj(t n ,Xj-i/2 -rf'(B))\ 
+ \uj(t„,x j - 1 / 2 - rf'(B)) - u e j_ 1 (t n ,Xj_ 1 /2 - rf'(B))\ 
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By Lemma 4.2 the first spatial derivative of |uf(i n )| is bounded by iV* j. Assume that \f"\ is bounded 
by S < oo. Since f'(B) < /3, for sufficiently small r, according to the inequalities (14 1 and (151, 



\u]{t n ,x 3 . 1/2 -Tf'(B)) - W |_ x (i„ ja;j _ 1/2 -rf{B))\ < De^h p+1+ ». 



Therefore, 



\A - B\ < Nlj5r\A -B\ + De ~» h p+1+t * . 
Due to the fact that r is very small, we have 

K(<-i/2) - ^0*7-1/2)1 =\A-B\<(1 + 2Nl J 8T)De p -'h p+1+ ». 
Now, plug ( |20[ ) into ( |19[ l, take the sum over all cells, we get 

(6,0n= (/(«*)- /(A &)n 

- B™ = ^ 1 [/(M p (a;7 +1/2 )) - /(it h (x7 +1/2 ))] C(^7+i/2) 

+ E™^ 1 [f{u{x+_ 1/2 )) - f(u h (xJ_ 1/2 ))\ £(x+_ 1/a ) 



+ 


(/O e ) - /K 










E717 1 


+1/2)) 


- /0A*7 4 


1/2))] f (^h 




/3 If U 2 (n) Ufa 


U 2 (o) 






+ 




1/2) 






+ 


/3E7T 1 | tt e (x+ 


-1/2) — 


U ( X j-l/2 


)l lf(4-i/a 


+ 


/3E7r 1 1 |it e (x7 


-1/2) — 


U i X j-l/2 


)l lf(*7-l/2 


+ 




-1/2) ~ 


h/ - 
U l - :r j-l/2 


)l If (4-1/2 


+ 




(fi) Ifa 






+ 


PY,™r 1 \u e (xj 


fl/2) - 


uP ( x j + l/2 


)l lf( S 7+l/2 



(22) 



By using ([22]), the estimates on the projection error given in Lemma 4.3 and the standard inverse 
inequalities ( |l0|. section 3.3), we can get computable constants C4,Cs and Ce,, such that 

d „ ^„ . C4 .,,11 , C5 , b+14-u , C*6 



f/f- 



llfllL 2( n) < ^||flk ( n) + ^ p+1+M 



h 



Integrating the last differential inequality, noticing the fact that £ = at r = 0, also noticing that 
r < r yh 1+a < jh 11 , we have a computable constant CV, such that 



llflk(n) < C 7 rh 



Since Hfllx,^) < v |fi| ||f ||l 2 (s2) , we have (16l and Lemma 4.5 proven. # 



Combining Lemma |4.4| Lemma |4.3| and Lemma |4.5[ we have the following theorem. 

Theorem 4.6 There is a computable constant J^S 1 ^), depending on the flux function f , the known 
constants of the interpolation/projection error estimates, the known constants of the inverse inequalities, 
and the components of the spatial smoothness indicator S„, such that 

\\u(t n+ i) - An+OlkCH) < r n h p+ ^(S V n)- 



4.2 Estimating u (t n +i) — u^ +1 , proof of dl2b 

The temporal smoothness indicator T% informs us about the boundedness of the temporal derivatives of 
u h at t — t n . We need to make sure that the boundedness of T„ can guarantee the boundedness of the 
temporal derivatives of u h for all t G [t„, t n +i]- 
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Lemma 4.7 There is a computable constant K, depending on the spatial smoothness indicator S p , such 
that 

ll« k (*n + r)|| too(n) < K|| W «) + Kh. (23) 
For each integer I 6 {1, • • • , k + 1}, there is a pair of computable constants ci and di, such that, for all 
t £ [t n , t n +i], 



Wggu \\ Loo ((i) < (l + cih a )\\— u h (t„)\\ Loo{(l) +d l h a . 
For each I, ci and di only depend on the -norms of the lower order derivatives. 
Proof. 

||w (t„ + T)\\ Loo (n) < \\u(t„ + r)|| Loo(f2) 

+ \\u(t„ + t) - u e (t n + r)|| ioo( n) 
+ \\u e (t n + t) - u p (t„ + T)|| Loo(fl ) 
+ \\u p (t„ + r) - u (t n + r)|| ioo (n). 



(24) 



II '"ll -too C^2) i s bounded by H'iiSillz.ro.), because the maximum of the 



4.5 



entropy solution u does not increase, 
and an application of the inverse 
The smallness of ||u — M e ||L 00 (n) 



4.3 



The smallness of \\u p — u Hi^m) can be obtained from Lemma 
inequality. The smallness of \\u e — u p || ioo (Q) is given by Lemma 

can be obtained by the same method used in proving (22 l. Consequently, we have the estimate (231 

In the proof of (24|, let's set the notation z = z(t) — u h (t„ + r), and z^ = j^jz. By differentiating 
the semi-discrete DG scheme 



0t,«)Of = {}{z),v x ) nj +f(z{x-_ 1/2 ))v(x^ 1/2 )~f(z(x :j+1/2 ))v{x j+1/2 ) 
with respect to t, we get 

(4 1) ,«)n 3 = {f'(z)z^,v x ) ai 

+ f'( z ( x J- 1/2 ))z {1) (xJ_ 1/2 )v(xf_ 1/2 ) 
- /' \z(xj +1/2 ))z w (xJ +1/2 )v(xJ +1/2 ), 



(4 2 \v)n 3 = (f'(z)z^,v x )n J +(f"(z)[z^] 2 ,v x )n J 

+ f' (z(xj_ 1/2 ))z {2) (xJ_ 1/2 )v(x+_ 1/2 ) + f" \z(xj_ 1/2 ))[z (1) (xJ_ 1/2 )] 2 v(x+_ 1/2 ) 
- f {z{x~ +1/2 ))z (2 \x~ +1/2 )v(x 7 +1/a ) - f"(z(xJ +1/2 ))[z (1) (xJ +1/2 )] 2 v{x~ +1/2 ), 

and similar equations for z">, I = 3, ■ • ■ , k + 1. It is easy to observe that the equation for z"' is linear on 
z^', I = 1,2, • • ■ , k + 1. Moreover, it depends on the derivatives of the flux function / and and products 
of lower order derivatives of z. 

In order to estimate the Loo-norm of z^\ we expand z^ 1 ' by the normalized Legendre polynomial 
basis functions {4>i,i : i = 0, • • • ,p} of Vh in each cell Slj. (4>j,i, (f>j,i)(ij = h/2. 

z (1) (r)| % =E? =o9 «(r)fe. 

Under this basis, let q' 1 ' be the vector consisting of all the q^} , j ' = 0, ■ • • , m — 1; z = 0, • ■ • , p. We can 
rewrite the equation for z^ 1 ' as 

hd (i) =4 (i) (r) (i) 

where A™(r) is the matrix obtained from the righthand side of the equation abo ut z ^ . The entries of 
j4^'(r) depend on the wave speed f'(z). Since ||j&||z, 0O (n) is bounded according to (231 and / is smooth, 
the entries of A^\t) are bounded. Besides, the entries does not depend on h, and there is at most 2p + 2 
entries in each row of A^(r) not equal to zero. Solving for from the last ODE, we get 

g W( r ) = e i/o TA(1) W* 9 W(0). 
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Figure 1: Smoothness Indicators, t=0.05 



From this solution, it is easy to see that, there is a constant A depending on the entries of A^(t), such 
that 

Ik^W-g^wlu < ^lk cl) (o)IU<M a lk cl) (o)l|oc. 

Due to the equivalence of \\z^' Hi^n) an d \\q ||oo> we have constants B = (p + \)^J (2p + l)/2 and 
C — V2, such that 

||* Cl) (T)|| Wn ) < \\z W (0)\\ Locm + \\z W (r)-z^(0)\\ Lxm 

< \\z a) (0)\\L 00 (n)+BA 7 h a \\qW(Q)\\ oc 

< ll^ (1) (0)|| ioc( n) + CBA 7 h a \\z^(0)\\L x( a) 



This proves (24 1 for 1 = 1, with ci = CBAj and di — 0. For I > 2, one can carry out the proof in the 
same way. # 



Remarks. (1) Lemma 4.7 confirms that one can essentially use the value of 4rrU h (0,t r , 



estimate of ^ J M ,l (r, t„, u%) in the entire time step [t n ,t n +i]. (2) The result of Lemma |4.7| serves our 



at 1 

purpose of local smoothness validation. However, neither the result nor the method of proof 3an/should 
be generalized to long term, because no numerical diffusion is taken into account. 



Based on the boundedness of the temporal derivatives proven in Lemma |4.7| it is trivial to conclude 
with the next theorem. 

Theorem 4.8 There is a computable function G{T^,S^), such that 

\\u h (t„+i) - u£+ilUi(n) < T* +1 g(T*, (25) 

5 Numerical evidences 



From the error estimation inequalities to ( 13 1, once we show the boundedness of all the components 



of the smoothness indicators, the rest of the error estimates is essentially a priori. Therefore, in order to 
demonstrate that our analysis works, it suffices to display the computed smoothness indicators. 

Example 1. In the first example, we solve Burgers' equation 

ut + {u/2) x = 
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Figure 3: Smoothness Indicators, t=2.0 



with the boundary condition uz,(t) = 1 and initial condition 

Ui(x) = 1 — (x/11) 3 sin(x) 

in x £ Q = [0, 10]. In this numerical example with p — 3 and k = 3, the cell size is ft = 0.05, while the 
time step size is r n = 0.005.. The solution has been computed in t £ [0,2]. The smoothness indicators 
at t = 0.05, £ = 1.05, and i = 2.0 are shown in Figure [T] Figure [2] and Figure [3] respectively. 

In each figure, the four plots in the top row are M^{— M n , M%, and M£, from left to right. The 
index j is dropped because each curve contains the values of M„j for j = 0, • • • , 199. The four plots in the 
second row are the temporal smoothness indicators Ut (t n ), Utt(t n ), utu(t n ), and u t ttt(tn)- The four plots 
in the third row are the jumps J°, J^, J^, and J^. In order to view the jumps from a better perspective, 
we show log h \J°\, log h log h and log h |J^| in the fourth row. Since J l nJ = D l nt] h p+1+tl - l(1+a \ 
the plot of log h \Jn j | = p + 1 + fj, — 1(1 + a) + log h |Z) l n A reveals the order (h ! ) of the jumps. Since p, fj, 
and a are all known, the values of D n j can be computed. Consequently, we can find D. 

It is easy to see the boundedness of M„ and T% in the figures when/where the solution is smooth. It 
is also easy to see that the order of the jumps J„ is as expected in the error analysis, or even smaller. 
These observations are sufficient to support the error estimates given in the paper. 

In addition, we have also observed some interesting phenomena. (1) logh\ J°\ — logn\ Jn\ ~ 2, logn\ Jn\ — 
log h \Jl\ w 1.4, log h \Jl\-log h \Jl\ w 1.8. There seems to be something related to the odd or even degrees 
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Figure 4: Smoothness Indicators, t„ = 0.0075, t = 0.12 



of the polynomials. (2) Long before the formation of a shock (t = 2.0, u x > —0.6), the fourth and 
third derivatives have grown significantly in a very narrow subdomain. The approximation benefit of 
the higher degree polynomials and the high order Runge-Kutta scheme will soon be lost locally at the 
spot. It seems that adaptive treatments need to kick in early. If not, there will be "numerical instability" 
showing up, ruining the numerical solution. 

In Figure[4] we show another numerical solution of the same problem, computed with h = 0.05 (same 
as before) and r n = 0.0075 (50% larger). The plots are made at t = 0.12, after 16 time steps from the 
initial time t = 0. With the improperly increased time step size, although the solution (presented by M° 
in the upper left corner plot) itself has not obviously shown anything wrong from the point of view of 
numerical stability (boundedness of solution, TVD, etc.), the higher order derivatives and jumps in the 
indicators have been increased significantly. The explanation is that the RKDG scheme for this problem 
with (p,k,h,r n ) = (3,3,0.05,0.0075) does not maintain numerical smoothness. As a consequence, the 
optimal approximation order must have been lost. The example seems to indicate the following: the 
strengthened CFL condition and the numerical diffusion from the Godunov flux are needed not only for 
numerical stability, but also for numerical smoothness maintenance. More attention should be paid to 
numerical smoothness when we are concerned with high order error estimates. 

The smoothness indicators can be used to diagnose the loss of numerical smoothness in an early stage, 
before too much damage is done to the global error. Of course, an algorithm needs to be designed for 
such diagnoses. We did run a separate case: after the first 5 steps at r n — 0.0075, T n is reduced back 
to 0.005. The spurious mode created in the first 5 steps were repaired in the following steps of smaller 
size. Nevertheless, the damage to the global error is done, unless we redo it. Further investigation in 
this direction can help in finding an optimal time step size. 

Example 2. In the second example, we show the solution of the Burgers' equation on [0, 10] with the 
initial condition 

Ui(x, 0) = - + -sin(7ra;/5) 

and the periodic boundary condition, k = 3, p = 4, h = 0.05, r n = 0.005. Figure [5] shows the numerical 
solution and its smoothness indicators at t = 1, when it is still far from any shock formation. 

In Figure [5] the five plots in the top row are M°(= u^), M%, and M*, from left to right. 
The five plots in the second row are the temporal smoothness indicators Un, Ut(t n ,), Ut t (t n ), ut tt (t n ), 
and Uf ttt (t n ). The five plots in the third row are the jumps J°, J*, J„ and J*. In the fourth row, 
we have log h | J° | , log h | J* | , log h \J„\, log h | J% | and log h | J* | . The values of the indicators are again 
what we expected and what we need to support the analysis. Obviously, the scheme has maintained the 
smoothness of the numerical solution, which guarantees that the local error of the next time step will be 
of optimal order. 



14 




Figure 5: Smoothness Indicators, p = 4, r n = 0.005, t = 1.0 



6 Conclusion remarks 

A. Choice of norm for error propagation analysis 

We prefer to use the Li-norm for error propagation analysis because of the well-known Li-contraction 
property. Other than the Li-contraction, a typical error propagation rate estimate for a time step con- 
tains a growth factor of the form 1 + CV. If we choose the Z/2-norm for error propagation, it is easy to 
show that the constant C is proportional to \J\u x f"{u)\. If we use numerical error propagation instead 
of PDE's error propagation, C will become bigger. "Bigger by how much" depends on the complexity 
of a numerical scheme. The appearance of u x f"(u) in the I/2-norm error propagation rate estimate 
implies that Z/2-norm error propagation analysis based on "worst case scenario" cannot be generalized 
to solutions with a shock or near a shock. Since large local error is expected to appear around the self- 
sharpening of a solution, the real scenario of a numerical solution is probably very close to the "worst 
case scenario". Li-norm error propagation analysis does not have this difficulty. 

B. How to deal with shocks and contact discontinuities? 

When there is a shock or contact discontinuity, it will be detected by the smoothness indicators, as 
shown in [3] . Certain quantitative criteria need to be developed to determine what kind of discontinuity 
is present according to the behavior of the indicators. It is also needed to determine if the discontinuity 
is well-caught, or some level of numerical "instability" has occurred. A decision should be made on the 
treatment of the discontinuity, including the use of a limiter or a local front tracking technique. After 
all of these have been done, we can consider error estimation. Error propagation is still to be estimated 
by using Li-contraction. Within each time step, in the smooth pieces of the solution, we can apply the 
error estimates given in this paper. At the discontinuities, we have to estimate the error according to 
the scheme. It is nice that the complexity of local error analysis does not get into the error propagation 
of the PDE. 

C. The process of sharpening before shock formation may be most difficult 

It might be the hardest to estimate error where a shock is forming but not yet fully developed. In 
this relatively wide space-time region, the solution's high order derivatives have become larger, causing 
difficulties for approximation. Adaptive algorithms need to be designed, and employed according to the 
smoothness indicators. As seen in Figure[3]of the first numerical example, the smoothness indicators can 
find the local sharp growth of the higher order derivatives and their jumps. The logarithm plots of the 
jumps have shown a clear exclusive pattern for a point of future shock. 
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D. Generalization to multi-dimensional problems 

We checked the proofs to the end of generalizing the results to 2-D scalar conservation laws. It seems 
to us that such a generalization should not meet any major difficulty. Generalization to hyperbolic sys- 
tems will face the lack of ii-contraction. 



E. a posteriori vs. a priori estimates 

The error analysis of this work is a posteriori because we depend on the computed smoothness indi- 
cators to compute the error estimates. However, if one can prove the boundedness of these smoothness 
indicators in advance, the error estimates can be converted to a priori error bounds. In this sense, under 
the concept of numerical smoothness, a prion and a posteriori error analysis has been united in the same 
framework. Moreover, our estimates are a posteriori in the sense that the smoothness indicators Sn and 

are computed after u c n has been obtained. As for the time step [t n ,t n +i], the smoothness indicators 
needed for the local error estimates of the step are computed before the local computation toward 
has started. In this sense, our error estimation is locally a priori, which will be more efficient if adaptive 
treatments are desired. 



F. Numerical smoothness of RKDG 

In the error analysis, we actually depend on the smoothness indicators to provide the needed numerical 
diffusion. That is, we take advantage of the RKDG method to include the needed numerical smoothness 
maintenance into the error analysis. The original designers of the scheme should get the credit for 
inventing a scheme with suc h pr operties. Since the numerical smoothness indicators S%, and T* are 



computed at (t„,Un), Lemma 
t G [t n , t„+i]. Lemma 



4.2 



4.2 



and Lemma 



4.7 



are needed to establish the smoothness of u and u h for 



shows the local smoothness preserving property of the PDE's strong solutions 
(in a special case useful for the analysis). Lemma 4.7 shows the local smoothness preserving property 
of the semi-discrete scheme. We only need these local smoothness proofs because smoothness is only 
needed in dealing with local error estimates. 
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